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Abstract 

o 

We consider the electromagnetic propagation in two-dimensional photonic crystals, formed by 
parallel dielectric cylinders embedded a uniform medium. The frequency band structure is com- 
puted using the standard plane-wave expansion method, and the corresponding eigne-modes are 
obtained subsequently. The optical flows of the eigen-modes are calculated by a direct computa- 
tion approach, and several averaging schemes of the energy current are discussed. The results are 
compared to those obtained by the usual approach that employs the group velocity calculation. We 
consider both the case in which the frequency lies within passing band and the situation in which 
the frequency is in the range of a partial bandgap. The agreements and discrepancies between 
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various averaging schemes and the group velocity approach are discussed in detail. The results 
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1/-) , indicate the group velocity can be obtained by appropriate averaging method. 
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I. INTRODUCTION 



When propagating through periodically structured media, i. e. photonic crystals (PCs), 
optical waves will be modulated by the periodicity. As a result, the dispersion of waves will 
no longer behave as that in a free space, and frequency band structures will appear. Under 
certain conditions, waves may be prohibited from propagation in certain or all directions, 
corresponding to partial and complete bandgaps respectively. The photonic crystals which 
could reveal bandgaps are called bandgap materials. 

Photonic crystals (PCs) are usually made of periodically structured materials which are 
sensitive to electromagnetic waves, and have been studied both intensively and extensively^, 
[J . The research of waves in periodic media has been first put on a theorematic basis 
stemmed from the efforts of Bloch[4] and Brillouinjs] towards electronic systems. An early 



survey was given in the excellent textbook[2, 6]. The general aspects of electromagnetic 
waves in photonic crystals are furthered reviewed in Refs. . A comprehensive survey of 
the literature can be referred to Ref. Q|. 

The propagation of electromagnetic waves in crystal structures is one of the top issues 
in the research of photonic crystals. The common theoretical approach to electromagnetic 
propagation in periodic media has been given in Ref. jjj, and may be summarized as follows. 
The Maxwell equations are first derived for waves in periodic media. By Bloch theorem, 
the solution can be expanded in terms of plane waves. The solution is then substituted into 
the governing equations to obtain an eigen-equation that determines the dispersion relations 
between the frequency and the wave vector that lies within the first Brillouin zone. These 
relations are termed as frequency band structures. Since it has been proved 3, |^ that the 
averaged energy velocity equals the group velocity which can be obtained as the gradient of 
the dispersion relations with the respect to the space of wave vectors, the investigation of 
electromagnetic propagation in periodic structures is thus reduced to the calculation of the 
group velocity from the band structures, thereby making a significant simplification of the 
problem. This approach may be termed as the group velocity approach and will be denoted 
as GVA hereafter. 

Since the average of the energy velocity, represented by the Poynting vector (S) ~ (E x 
H), is performed over the whole unit cell of periodic media, a few questions may be asked 
about this group velocity approach. The first question is whether such averaged energy flow 
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can fully depict the actual electromagnetic energy flow in periodic media. Second, since in 
actual measurements it is often hard to detect the physical quantities within the periodic 
media, how to obtain the group velocity without having to put a detector into the media 
needs also to be considered. Deducing the group velocity from practical measurements is in 
fact an important task in the photonic crystal research. More explicitly, the question is how 
to relate practical measurements with the group velocity which can be derived theoretically 
from the band structure calculation. 

The purpose of the present paper is to examine the question of how well GVA can describe 
the actual electromagnetic (EM) energy flow in periodic structures, and discuss the issue 
how to obtain the group velocity from practical perspectives. We will compare the averaged 
energy flow obtained from GVA with the results from the direct computation of the energy 
current. Various schemes are proposed to obtain the averaged current and compared, so to 
find the most appropriate schemes in deducing the group velocity in practical measurements. 

To simplify our discussion yet without losing generality, we will only consider the prop- 
agation of electromagnetic waves in two dimensional periodic media, i. e. two dimensional 
photonic crystals, which are made of arrays of dielectric cylinders. A reason of using this 
type of systems is that the systems are experimentally ready, and therefore the conclusion 
derived from the present paper can be verified. The significance of the present research 
is two folds. First, it provides useful information about how to obtain the group velocity. 
Second, it shows under what conditions the group velocity approach is valid. 

The paper is structured as follows. In the next section, we will present the usual formu- 
lation of Maxwell's equations for EM waves in periodic structures. The energy flow will be 
formulated, and the GVA will be outlined in general. In section III, the particular system 
will be discussed, and a few averaging schemes of the energy current will be proposed. The 
numerical results for a number of situations will be presented in Section IV. The paper will 
be concluded by a summary in section V. The proof of the equivalence between the aver- 
aged energy velocity and the group velocity will be briefly repeated in the Appendix, for the 
readers' convenience. 
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II. THE GENERAL FORMULATION 



A. Bloch wave solution, energy current, and energy density 

The EM waves in two dimensional media can be separated into two cases of polarization: 
(1) the s polarization or the E-polarization, that is, the electric field is along the z axis 
perpendicular to the plane of propagation, which is the plane of the periodicity; and (2) 
the p polarization or the H polarization, with the magnetic field being perpendicular to the 
plane of propagation. Here we outline the determining equations for EM propagation in two 
dimensional periodic media. The details can be referred to in Refs. 

For both polarizations, the governing equation can be unified as 

V Qvp(rO) + -^Vif) = 0, (1) 

where c 2 , p and q are two dimensional periodic functions, depending on the properties of 
the medium. For the s polarization, p stands for E z , with p(r) = p(r) and q = p(r). For 
the p polarization, p denotes the magnetic field H z perpendicular to the wave propagation. 
In this case, p(r) = e(r) and q = e(r). 

Due to the periodicity, we can make the following expansions 

- = °{G)e id - ? , and ±- = £ x(G)e iGV , (2) 
p / — ' c z q ' 

G G 

where G are the reciprocal vectors. 

By Bloch's theorem, the solution to Eq. can be expressed as 

p R (r) = e^uz(r), (3) 

where K is the Bloch vector that lies within the first Brillouin zone, and is a periodic 
function with the periodicity of the medium; therefore pg is the eigen field corresponding 
to the Block vector K. The function can be expanded as 

u R (f) = J2 A K(G)e^. (4) 

G 

Substituting Eqs. (J2J) and (jSJ) into Eq. (JTJ), we obtain 

[<r(Gi) ((K + G 2 ) -(K + G 2 + &)) - u; 2 x(Gi)] A R {G 2 ) = (5) 

Gl 
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From this equation, we can find a secular equation that determines the dispersion relation 
between ui and K, 



det ({K + G 2 ) ■ {K + G 2 + G,)) - uj^gA _ = 0. (6) 

V /J Gi,G2 

Once the dispersion relation is determined, and the coefficients Ag(G) can be obtained from 
Eq. (jSJ). The EM waves can be subsequently obtained from Eqs. (j2J) and (J1J). 

When either the electrical E or magnetic field if is determined, corresponding to the s 
or p polarization respectively, the magnetic or electric field for the Bloch vector K can be 
determined from 

V x = —iuoeEg, or V x Eg- = iuj/jHg, (7) 

where we have assume e~ luJt time dependence. 

By Eq. 0, the time averaged flux of energy at any spatial point is subsequently obtained 
from 

4 = 1 -MEk x H*g\, (8) 
where refers to the complex conjugate operation. And the time averaged energy density 
is 

U li = 1 -[e\E l} \ 2 + f ,\H l} \ 2 }. (9) 
From these two equations, we can directly calculate the EM energy flow and density. 



B. Group velocity approach (GVA) 

In principle, the EM propagation in periodic media can be inferred from a direct calcu- 
lation using Eqs. (jSJ) and Q. However, it is common to use the group velocity method to 
discern the EM energy flows in periodic media. This approach is based on the following 
theorem^]. The group velocity in periodic media is defined as 

v 9 = V^u;. (10) 

Or in another form, 

5oo = v g ■ 5K. (11) 
The energy velocity in the periodic media is defined as 

( e) " ±fud* r -( Ujt y (12) 



where V is the volume of a unit cell, and the integration is performed in the unit cell. In 
two dimensions, this equation is reduced to 

{e) ~ilu g d'r-{u g y (13) 

with S being the area of a unit cell and the integration being restricted to the unit cell. It 

rj 

can be proved [2] that 

(v e ) =v g = V R u. (14) 

This equation will be referred to as the equivalency theorem. 

By Eq. (JHJ), the task of finding the EM propagation in periodic media is reduced to 
calculating the group velocity, i. e. v g = V gu. This procedure may be called the group 
velocity approach (GVA). As the group velocity is relatively easy to be calculated from the 
band structures, this method has been widely used. 

From Eq. (|13jl. we see that the spatial average of the current is taken over the whole 
area of the unit cell for two dimensional cases. Experimentally, the energy flux is normally 
determined by measuring currents flowing through certain surfaces. In the following, we will 
examine whether and when the spatially averaged currents can represent the actual flows. 
We will examine a few averaging schemes. 



III. THE SYSTEMS AND THE VARIOUS AVERAGING SCHEMES 

In the following, we will compare the energy flux obtained from the GVA with the results 
from the direct computation of the current given by Eq. ©• 



A. The system 

The systems considered here are two dimensional photonic crystals made of arrays of 
parallel dielectric cylinders placed in a uniform medium, which we assume to be air. Such 
systems are common in both theoretical simulations or experimental measurements of two 
dimensional PCsQ. For brevity, we only consider the E-polarized waves (TM mode), that 
is, the electric field is kept parallel to the cylinders. The following parameters are used in the 
simulation. (1) The dielectric constant of the cylinders is 14, and the cylinders are arranged 
to form a square lattice. (2) The lattice constant is a and the radius of the cylinders is 
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0.3a; in the computation, all lengths are scaled by the lattice constant. (3) The unit for 
the angular frequency is l-ncja. After scaling, the systems become dimensionless; thus the 
features discussed here would be applicable to a wider range of situations. 

For the E-polarization , i. e. the electrical field points along the z axis, the axis of the 
cylinders, Maxwell's equation is simplified as 

[V 2 + e(fV 2 /c 2 ] E(f) = 0, (15) 

where the electrical field is represented as 

E(f) = e- iut E{f)z. (16) 

By Bloch's theorem, we have 

E(r)=e^J2 E d^^ (17) 

G 

where the wave vector K is in the first Brillouin zone, and G is a reciprocal vector. Substi- 
tuting Eq. (JT7j) into (JT3J), we obtain the usual matrix equation 

~oo 2 /c 2 -(G + K) 2 

G> 

where 

B(G - G') = 27r(e - l)aJi(|G - G'\)/{\G - G'\d 2 ), 

with Ji being the Bessel function of the first order. After E is obtained, the magnetic field 
can be deduced by employing Eq. ((7j). The energy current can be therefore derived from 
Eq. ©. 

B. Various averaging schemes 

To compare the energy current determined from the GVA with that obtained from the 
direct computation method, we take the following procedure. Due to the periodicity and 
symmetry, it is sufficient to just consider the energy current in a unit cell, which takes the 
square shape. For a given frequency, a Block wave vector K can be determined from the 
secular equation for the band structure in Eq. (fTSj) . The group velocity, therefore the energy 
velocity, will then be calculated for K from the band structure as 

Vg = VgU. (19) 
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- ^2 u?B{G — G')Eq,, 



(18) 



This consideration is illustrated by Fig. ^ Here, the coordinates are shown in the figure. 

According to the aforementioned equivalency theorem, the direction of the group velocity 
will be the direction of the spatially averaged EM energy current, i. e. 



(J) = 1 / d 2 rJ. (20) 

J S 

We note that the integration is performed over the whole area of the unit cell, which includes 
the areas occupied by the scatterers. In actual experiments, it is often difficult to probe 
the currents within the areas taken by the scatterers, therefore we may replace the whole 
integration by a partial integration 

(J') = i / d 2 rf. (21) 

o S 1 

Here the integration is performed over the area 5" which excludes the areas occupied by the 
scatterers. 

In practice, the EM energy current through this unit cell can also be calculated from the 
flux across the two sides of the cell denoted by AB and BC, i. e. 

I y ~ / dxJ • y, and I x ~ / dyj ■ x. (22) 
Jab Jbc 

Here J is given by Eq. (JHj). 

We will compare the results in Eqs PU|). (j2*Tj) and (|2*2*|) with that obtained in Eq. (JT4*j) by 

the GVA: 

Iy VA ~ (Vj^) ■ y, and ~ (V^) • x. (23) 

In the present paper, we label the averaged current in Eq. as case 1, that in Eq. (J2UJ) 
as case 2, and that in Eq. (|2*T|) as case 3. 

To simplify our discussion, we will compare the ratio between the averaged current in 
two directions. 

• The GVA case. The angle of the group velocity is determined as 

^ = tan-f^V (24) 



Case 1: y- from Eq. (1221) . we represent the ratio by the angle 



h = tan- 1 rf ) , (25) 



S 



{J)y 

" 1IU111 J-l/4. Il-A ^) . 111C UUIICDJJU1 

-1 / (J) 



Case 2: y-=p from Eq. (|20j) . The corresponding angle is 



Case 3: from Eq. (|21jl . The corresponding angle is 

(J )i 



11 Will J_/^. II- ± I) ■ 111C ^UllCDj^ 

.-. ( ( J> ) 



There are other options in obtain the averaged current. We refer to the setup shown in 
Fig. If the detection is along line AB, the averaged current vector may be obtained as 

1 r B - 

{Jab) = ~r- / dxJ. (28) 

J-'AB J A 

If the detection is on line BC, the averaged current vector will be 

1_ rB 

^CB Jc 

Here Lab and Lqb denote the lengths of the two sides of the unit cell. We may also 
consider the sum of these two averaged current if the detection is made on both AB and 
BC. Correspondingly, there are other three possibilities 

• Case 4: ffsh. from Eq. (j2~9l . The associated angle is 

(JcB)x 

04 = tan -i (i^k I . f30) 



(Job) = ~r- I dyJ. (29) 



Case 5: ^sk. f rom Eq. J29. The associated angle is 

(JAB)* 



1IU111 J-J4. ■ 111C CIDOU^J 

-1 / (•/as) 



5 = tan- 1 ff^ ] . (31) 
\{Jab)x t 

Case 6: ^^±^4r from Eqs. (El and (El. The associated angle is 

{{Jab)+{Jcb))x 



tan- f «fo + ^ 



(32) 



Since the energy or the current fields in the areas occupied by the scatterers, may not be 
easy to detect, thus there is another possibility in case 4, 5 and 6. That is, the contribution 
from these areas are excluded. Later we will compare the results obtained from various 
averageing schemes. 



IV. THE RESULTS AND DISCUSSION 



The frequency band structure is plotted in Fig. |21 A complete band gap is shown between 
frequencies of 0.22 and 0.28. Just below the complete gap, there is a regime of partial band 
gap in which waves are not allowed to travel along the TX or [10] direction. We will consider 
waves two frequencies: one is at 0.16 which is in the first passing band, and the other is at 
0.19 which is within the partial bandgap. 

A. Two dimensional imaging of energy and energy current fields 

First we study the spatial behavior of the energy density fields and the local flows of the 
eigen-modes which are characterized by the Bloch wave vectors. The current is computed by 
Eq. (JBJ), while the density field is calculated by Eq. (JUJ). The results are shown in Figure El 

The left and right panel of Fig. El describes the result for frequency 0.16 and 0.19 re- 
spectively. The former lies within the first passing band, whereas 0.19 is within the first 
partial bandgap regime. Within this partial band, the waves are forbidden from propagation 
along the TX direction, i.e. [10]. For each frequency, we have considered three eigen modes 
represented by three Bloch vectors which are given in the figure caption. The two principle 
directions of the unit cell, TX and FM, are also shown in the figure. 

Here we observe the following. First we discuss the case of the passing band in (al), (a2) 
and (a3). Overall speaking, the flows of the energy indicated by the black arrows tend to 
flow along the direction indicated by the Bloch vectors. This feature is more obvious for 
the local current within the dielectric cylinders. For the frequency within the partial gap, 
however, we observe that most of the current flows may not be along the direction of the 
Bloch vector. Fig. bl) shows that for small angles with reference to the [10] direction, all 
the local currents tends to flow along the direction which is close to the direction of [01]. 
That is, the energy flows are nearly vertical, but they cannot be exactly vertical, as the 
direction of [01] is a forbidden direction. 

When the direction of the Bloch vector is tilted more and more away from the direction 
[10], an interesting feature prevails. That is, the currents eventually tends to flow into 
the direction of TM, i. e. the [11] direction. This feature is clearly demonstrated by the 
examples in Fig. 3(b2) and (b3) for which the Bloch vector points to the angles of 30° and 
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45° respectively. 

Comparing the results for the passing band in the left panel and the results for the 
partial band gap in right panel of Fig. El we may conclude that the electromagnetic flows in 
periodic structures or photonic crystals will highly depend on the band structures. There are 
significant differences in the current behavior between the situation in which the frequency 
is located in a passing band and the case in which the frequency is within a partial band gap. 
The result shown by Fig. Efb2) suggests that an effect of the partial band gap is to bend the 
current towards a direction which is to avoid the forbidden direction as much as possible; 
in the present case, it is the direction of TM or [11]. Such feature may render possible new 
applications of partial band gaps in manipulating EM waves in opto-electronic devices. The 
results in Fig. |3] also shows that both the energy fields and the current fields are not uniform 
inside the unit cell. This feature indicates that the energy velocity is not also not uniform. 
The energy is more concentrated within the regimes occupied by the scatterers. 

B. Comparison of different methods obtaining the averaged currents 

The results in Fig. |3] indicate that the local energy current is not uniform within a unit 
cell of a periodic structure. How to describe the overall energy flow in such a non-uniform 
situation of periodic structures thus poses an important issue. 

As described in Section lHB| we mentioned that the common theoretical approach is based 
on the equivalence theorem between the group velocity and the averaged energy velocity^]. 
The average is taken over the whole volume in three dimensions or the whole area in two 
dimensions, which include the volumes or the areas occupied by the scatterers. In actual 
experiments or observations, however, it may be difficult to probe the whole volume or the 
whole area to deduce the information of the averaged energy current. In particular, the 
currents or density within the volumes or the areas occupied by the scatterers are hard to 
detect. As matter of fact, there is no report on detecting the energy or energy current over 
the whole volume or the area to our best knowledge. The usual detection is made either at 
one particular spatial point or on a certain surface. In addition, it is often the intensity field 
that is measured. 

In this section, we will compare the results obtained from the various averaging schemes 
outlined in Section IIIIBI The results will be compared with those obtained by the GVA. 
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For brevity yet without losing generality, we consider the two frequencies from Fig. 0.16 
and 0.19; one is in the passing band and the other one is within the partial bandgap regime. 

First, we compare the first three cases of averaging and the GVA. The three scenarios 
are given by Eqs. (J25|) . (|26|) . and (J27|) respectively. The results are shown in Fig. 0] Here 
the directions of various averaged currents and the group velocity are plotted against the 
direction of the Bloch wave vectors. The directions are represented by the angles of the 
corresponding vectors with reference to the x axis. 

For either the passing band or the partial bandgap cases, we see that the results from 
the averaging scheme 2, that is the average is taken over the whole area of the unit cell, 
fully agree with the results from the GVA. This verifies the equivalence theorem. It can be 
also seen that the results from the scheme 3 also agree reasonably well with the GVA. This 
implies that as long as the current inside a periodic medium can be measured, the group 
velocity can be well deduced by the averaging scheme, no matter whether the areas occupied 
by the scatterers are excluded or not. The situation is noticeably different for the averaging 
scheme 1, in which the information relies on the measurements along the two boundaries of 
the unit cell. This scheme can reproduce the results of the GVA considerably well for the 
passing band case in Fig. EJa), but fails for most of the Bloch wave vectors in the partial 
bandgap case in Fig. |3Jb). In the later case, the agreement recovers as the angle approaches 
45 degree, i. e. as the direction of the Bloch vector approach that of TM. 

We have also compared the three other averaging schemes. The results are presented 
in Fig. 0D Here, we have considered both the situation in which the areas occupied by the 
scatterers (cylinders) are included and the case in which these areas are excluded. Here we 
observe the following. (1) For both the passing band and partial bandgap case, the average 
over any single side of the photonic crystal will either overestimate (case 4) or underestimate 
(case 5) the angle of the group velocity, no matter whether the areas of the cylinders are 
included or not. In other words, the group velocity will not represent the energy current 
averaged only along one observation line, as represented by the results in case 4 and 5. 
Relatively speaking, when the contributions from the areas of cylinders are included, the 
results move closer to that obtained from the GVA. (2) The results from scheme 6 are in 
excellent agreement with the results from the GVA, no matter whether the areas of the 
cylinders are included or not. This feature implies that this averaging scheme is a good 
candidate in inferring the group velocity of photonic crystals when the energy density and 
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energy current inside the crystals cannot be readily probed. 

Some common features can be discerned from Figs. 0] and 03 In the passing band case, 
the direction of the group velocity is close to that of the Bloch vector. The relation between 
the angles of the two is almost linear, referring to Fig. E(a) and Fig.|5Jal) and (bl). For the 
first partial bandgap situation, the angle of the group velocity decreases as the angle of the 
Bloch vector K increases. Actually as long as the angle of the Bloch eave vector exceeds 25 
degree, the angle of the group velocity saturates to the value of 45 degree. This indicates 
that partial gaps tend to bend currents into certain directions, allowing for possible novel 
applications of photonic crystals in the partial bandgap regimes joj]. We have done further 
simulations, and found that these features are also true for other frequencies in the first 
passing band and first partial bandgap. 

C. A note on the energy density and energy flow in photonic crystals 

From Fig. 01 we see that the spatial distribution of the energy density and the energy 
current is highly non-uniform in the unit cell. This will also indicate that the local energy 
velocity, defined as ? e e i with Jg and Ug being given in Eqs. (JBJ) and ©, is also be 

K 

non-uniform. This will have some implications on the observation of the energy density and 
energy flows. Since the current J equals Uv e , then with a given current magnitude, a larger 
local energy density (intensity) implies a smaller velocity. For instance, consider two current 
vectors, J\ and J 2 , with the same magnitude, but perpendicular to each other. Clearly, the 
summation of the two vectors, Jt = J\ + J2, will give a total vector which lies between the 
two vectors. If the magnitudes of the two corresponding current velocities are not equal, 
then the larger is the velocity magnitude, the smaller will be the energy density. As a result, 
the apparent energy density field will not be aligned along the direction of the total current 
Jt- This implies that the current flow may not be readily obtainable by just measuring 
the energy intensity field. We may also look at this problem from another perspective. 
The current flow deduced from the group velocity approach, or the spatial averaged current 
method in case 2 and 3 say, may not be able to describe the apparent energy intensity field 
which is actually the quantity to be measured. 
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V. SUMMARY 



We have considered the electromagnetic propagation in two-dimensional periodic arrays 
of dielectric cylinders embedded a uniform medium. The frequency band structure is com- 
puted using the standard plane-wave expansion method, and the corresponding eigne-modes 
are obtained subsequently. The spatially dependent optical flows of the eigen-modes are cal- 
culated by a direct computation approach. A few averaging schemes for the energy flows 
are discussed. The results are compared to those obtained by the common group velocity 
approach (GVA) which is based upon the group velocity calculation. We have considered 
both the case in which the frequency lies within passing band and the situation in which 
the frequency is in the range of a partial bandgap. It is shown that some average schemes 
may reproduce well the results of the GVA. With these schemes, the group velocity can be 
deduced in measurements. The research provides a useful information about how to obtain 
the group velocity and what information the traditional GVA can provide. 
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APPENDIX A: PROOF OF v e = v g 

For the convenience of the reader, we present a proof of v e = v g . 
In photonic crystals, the waves assume the Bloch form, 

E(r) = Eg(r)e^, (, 

and 

H(f) = Hg(f)e^, 

where Eg and Hg are periodic functions, and K is the Bloch wave vector. 
By substituting Eqs. (|A1|) and (|A2|) . we have 

^ 2Re[(Eg x Ht}} 

v e = =; — , U 

(^| 2 + /^| 2 > 

where (■) denote the integration over the unit cell. 

Taking eqs. (jAlj) and ()A2j) into the Maxwell equations ((Zj), we have 

V x Hg + iK x i?^ = —iueEg, (i 

and 

V x + iK x Eg = iufiHg. (i 

Suppose now that K is changed by an infinitesimal amount, SK. This will induce chan 
5u, 5 Eg, and SHg, and they are related as 

V x 5Hg + i8K x Hg + iK x <5i?^ = —ibuotEg — iue5Eg (i 

and 

V x + i5K x £^ + iK x 5Eg = iSoofiHg + iufiSHg. (j 

In the following, let us ignore the subscript K. 
Then we have 

VxH + iKxH = -iueE, (j 

and 

VxE + iKxE = iuj^iH, (j 

and 

V x 5H + i5K x H + iK x 5H = -i5iueE - iuje5E, (A 
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and 

V x 5E + i5K x E + iK x 5E = i5ufiH + iuj^5H. (All) 

Also we have 

V x H* -iK x H* = iueE*, (A12) 

and 

V x E* -iK x E* = -iufiH*. (A13) 
Multiplying Eq. (|A"To|) with E*, we get 

E* ■ (V x SH) + iE* ■ (5K x H) + iE* ■ (K x SH) = -i5uje\E\ 2 - iueE* ■ SE, (A14) 
and its complex conjugate is 

E - (V x «T) + z£ ■ (5/? x /T) +iE-(K x = -«We|£| 2 - ztueE • (A15) 
Similarly, multiplying Eq. ()A11|) with H*, we have 

H* ■ (V x 5£) + z'F* ■ (5K xE) + iH* ■ (K x SE) = %5u^\H\ 2 + iu/iH* ■ SH, (A16) 
and its complex conjugate 

# • (V x + iff ■ x £*) + iH -(K x SE*) = iSujfi\H\ 2 + iu;/ff? ■ 5ff*, (A17) 
Subtracting Eq. (jA14|) from (jA16|) and using 

A- {B x C) = B ■ (C ■ A) = C ■ [Ax B), (A18) 

we get 

iSoo(e\E\ 2 + /i|ff | 2 ) + iu(eE* ■ SE + /iH* ■ SH) = iSK ■ (E x H* + E* x H) 
+H* -[iK x SE + V x SE) - E* ■ [ik x SH + V x SH] . (A19) 

The last a few terms on the RHS of Eq. (jA19|) can be derived as 
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H* ■ [V x 5E + iK x SE] — E* ■ [iK x SH + V x 5H] = 
+5H ■ (iK x E*) -l*-(Vx <5#). 



-522 • (itf x #*) + F-(Vx &E) 

(A20) 



Now we use 

iK x #* = V x H* - iueE*, (A21) 

and 

iK x E* = V x E* + iufiH*. (A22) 

Then 

I = -8E ■ (iK x H*) + H* ■ (V x + 5# ■ (iK x E*) - E* ■ (V x 
= Kje£* • - 5£-(Vx #*) + H* ■ (V x &E) 
+ iu;//# ■ «T + 5H-(V x E*) - E* ■ (V x <5#) 

or 

/ = io;e£* • 5£ + iu/iH -5H* + V ■ (#* x 5#) + V ■ (8H x E*). (A23) 
Therefore Eq. (|A"T9j) becomes 



«W(e|£| 2 + fi\H\ 2 ) + zo;(e£* • SE + fiH* ■ 5H) 



iSK ■ (E x H* + E* x H) + 

iujeE* ■ 5E + iufiH ■ 5H* + 

V • (H* x 8E) + V- (5H x E*). 

(A24) 



Its complex conjugate is 



- i5uj(e\E\ 2 + fi\H\ 2 ) - iu(eE ■ 5E* + fiH ■ 5H*) = -i6K ■ (E* x H + E x H*) + 

-iuotE ■ 5E* - iufiH* ■ 5H + 
V- (H x 5E*) + V ■ (5H* x E). 

(A25) 



17 



Due to the periodicity, for any periodic function A we have 

J dvV -A = J dS-A = 0. (A26) 

Now we subtract Eq. (|A25J) from Eq. (|A24J) . then we perform the volume integration over 
a unit cell to get 

2i5uo((e\E\ 2 + fi\H\ 2 )} = AidK ■ Re[(E x #*)]. (A27) 

That is 

5cu = v e ■ 8K. (A28) 
Comparison between Eqs. (fTT|) and (|A28|) leads to 

v g = v e . (A29) 

From the above derivation, we understand the conditions for v g = v e to be held are that 
(1) the variation 8K can be arbitrary; (2) when the variation 5K is small, the changes in u, 
Ek and Hk are also small. 
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X M 

Band Structure 



FIG. 2: The band structure of a square lattice of dielectric cylinders. The lattice constant is a and 
the radius of the cylinders is 0.3a. TM and YX denote the [11] and [10] directions respectively. 
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FIG. 3: The imaging of the intensity fields \E\ and the current of eigen- modes. Two frequencies 

are taken: 0.16 and 0.19 for the left and right panel respectively, (al) K = (0.77ir/a, 0), i.e. the 

Bloch vector points to the angle of 0°; (a2) K = (0.667r/a, 0.387r/a), i.e. the Bloch vector points 

to the angle of 30°; (a3) K = (0.547r/a, 0.547r/a), i.e. the Bloch vector points to the angle of 

45°. (bl) K = (0.997r/a,0.42vr/a), i.e. the Bloch vector points to the angle of 23°; (b2) K = 

(0.877r/a, 0.5l7r/a), i.e. the Bloch vector points to the angle of 30°; (b3) K = (0.697r/a, 0.697r/a), 

i.e. the Bloch vector points to the angle of 45°. 
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(a) Passing band 
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(b) Partial band gap 
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FIG. 4: Comparison of four different ways obtaining the averaged EM current in a unit cell. The 
three cases are from Eqs. (|25j). (|2f)|) and (|2Tj) respectively. The results from the GVA are obtained 
from Eq. (T2IJl . 
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(a1) Passing band 
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(a2) Partial band gap 
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(b2) Partial band gap 
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FIG. 5: Comparison of the results from the averaging schemes in case 4, 5 and 6. The results in the 
left panel are obtained as the areas occupied by the cylinders are excluded in the average, whereas 
the results in the right panel are obtained when the areas are included. The three schemes are 
from Eqs. (|31|). (|3U|) and ()32l) respectively. The results from the GVA are obtained from Eq. (|24|) . 
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